Thermodynamics of the dissipative two-state system: a Bethe Ansatz study 
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The thermodynamics of the dissipative two-state system is calculated exactly for all temperatures 
and level asymmetries for the case of Ohmic dissipation. We exploit the equivalence of the two- 
state system to the anisotropic Kondo model and extract the thermodynamics of the former by 
solving the thermodynamic Bethe Ansatz equations of the latter. The universal scaling functions 
for the specific heat C a (T) and static dielectric susceptibility \a (T) are extracted for all dissipation 
strengths < a < 1 for both symmetric and asymmetric two-state systems. The logarithmic 
corrections to these quantities at high temperatures are found in the Kondo limit a — *• 1 _ , whereas 
for a < 1 we find the expected power law temperature dependences with the powers being functions 
of the dissipative coupling a. The low temperature behaviour is always that of a Fermi liquid. 
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I. INTRODUCTION 

The low temperature properties of many physical sys- 
tems, such, as quantum tunneling between flux states in 
a SQUID, El, two- level atoms coupled to the electromag- 
netic field in quantum optics,cl, or tunneling dislocations 
and point defects in solidsd can bft-described by the dis- 
sipative two-state system modelcrO. All these systems 
have the common feature that at low energies the sub- 
system investigated (the flux states of the SQUID, the 
atom, or the crystal defect in the examples above) can 
occupy only two distinct quantum states, which are cou- 
pled to a continuum of excitations (the electromagnetic 
field, the phonons or the conduction electrons^, leading 
to the appearance of dissipation in the systems. 

In the dissipative two state system (DTSS) model the 
two distinct quantum states above are described in terms 
of a pseudospin <7j, with a z = ±1 corresponding to the 
two states of the subsystem. Generally, these two states 
have slightly different energy with an energy difference e 
(also called asymmetry energy) and the decoupled two- 
state system (TSS) can tunnel between them with a tun- 
neling amplitude A. The heat bath is modeled by a 
continuum of independent quantum oscillators with den- 
sity g(u>) coupled linearly to a z with a frequency de- 
pendent coupling g(u)). A detailed analysis shows that 
the dynamical properties of the two-state system are 
uniquely determined by the environment's spectral func- 
tion, J(oj) ~ g(u>)g 2 (iv) ~ lu s U. 

Here, we only concentrate on Ohmic dissipation cor- 
responding to s = 1, i.e. J (to) ~ w. This includes, 
for instance, the important case of a tunneling defect 
in a metal, where the low-energy bosonic excitations are 
electron-hole pairs close to the Fermi surface. These exci- 



tations have a linear dispersion, and in a first approxima- 
tion their coupling to the defect is energy-independent, 
leading to J(oj) — 27raw for to < uj c with uj c a high-energy 
cutoff in the model of the order of the Fermi energy Ep. 
The dynamical behaviour of the TSS as a function of the 
coupling strength a has been the subject of extensive 
studies during the past two decades and it is well under- 
stood by nowu. To distinguish between the different cases 
it is useful to introduce the zero temperature spin correla- 
tion function S(t) = $t{a z (0)a z (t)) . (a) For a < 1/2 the 
TSS oscillates between the two states a z = ±1 and S(t) 
has an oscillatory behaviour. However, the environment 
introduces some decoherence in the system, reflected in 
the exponential decay of the envelope of S(t). It also 
renormalizes the tunneling amplitude: A — > A. r < A. 

(b) In the parameter range 1/2 < a < 1 the coher- 
ent oscillationsa become completely suppressed, and S(t) 
shows an exponential behaviour without a change of sign. 

(c) Finally, for a > 1 and a finite level asymmetry the 
TSS becomes localized in the lowest quantum state (at 
T = 0)lij. In this case S(t) tends to a finite value as 
t — > oo. It is important to note that the localized state 
obtained is immediately destroyed once assisted tunnel- 
ing oij-assistcd pair tunneling is included in the DTSS 
modeOB 

In the present paper we study the thermodynamics of a 
dissipative TSS model in the parameter range < a < 1. 
To this purpose we exploit a mapping between the spin 
anisotropic Kondo model (AKM) describing an impu- 
rity spin coupled to the spin density of the conduction 
electrons via an anisotropic exchange interaction (see 
Sec. lib) and the dissipative-JSS modelpcl, and study 
the Bethe Ansatz equationstU for the former both nu- 
merically and analytically. As discussed in Sec. [IB and 
Appendix |A], within this mapping the tunneling ampli- 
tude maps to the spin flip scattering amplitude, A «-> Jj_, 
the asymmetry energy e corresponds to a local magnetic 
field h applied to the impurity spin, and the dissipation 
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strength a is related to the coupling J z in the Kondo 
model. It is very remarkable that the a values separat- 
ing the three different regions of the DTSS model are 
mapped to some special points in the parameter space 
of the Kondo model. The point a = 1 turns out to 
correspond to the case J z — 0, separating the ferro- 
magnetic (J z < O a > 1) and the antifcrromagnetic 
(J z > <^> a < 1) regimes in the Kondo model. While in 
the first case the Kondo model scales to a finite fixed 
point, in the second the Kondo fixed point turns out 
to be at infinite coupling. The crossover to the 'strong 
coupling' regime happens at the so-called Kondo energy, 
T K JiSEl which can be identified with the renormalized 
tunneling amplitude in the DTSS model, A r ~ Tr. The 
other special point, a = 1/2, can be shown to be equiv- 
alent to the Toulouse linetL!3 of the AKM. Along this line 
the Bethe Ansatz (BA) equations simplify enormously, 
and the model can be described by a simple resonant 
level model without interaction. 
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FIG. 1. (a) The specific heat C(T) on a linear tempera- 
ture scale for the symmetric case at weak and strong dissi- 
pations. The position of the maximum of the Schottky peak 
in C(T) is of order A r but its exact location changes with 
a, as does the peak height, (b) The quantity A r C(T) /ksT 
shows non-monotonic behaviour for weak dissipations, and 
monotonic behaviour for a > 1/3. 



<x=l/5 
<x=4/5 




We give a detailed description of the thermodynamics 
of the DTSS model, solving the BA equations for ar- 
bitrary asymmetry and temperature for both a > 1/2 
and a < 1/2. For the sake of completeness and clarity, 
we also included a detailed analysis of the two models 
to demonstrate, how the different concepts such as scal- 
ing, strong coupling limit, energy scales, etc. appear in 
the AKM and the DTSS model. As we shall see, there 
are clear indicators in the specific heat, C(T), for dis- 
tinguishing weak from strong dissipation limits. This is 
not evident in C(T) directly, which shows a Schottky 
anomaly at fc^T ~ A r for all dissipation strengths a < 1 
as depicted in Fig.|l|a . However, as seen in Fig. |b, and 
as we shall discuss in detail later, the quantity C(T)/T, 
shows quite different behaviour at weak and strong dis- 
sipations. For dissipations a < 1/3, with no asymmetry, 
C(T)/T is found to have a peak at ksT ~ A r indicating 
the expected tendency towards activated behaviour of the 
two-level system as a — > 0. A quite different behaviour is 
found for a > 1/3 where we find that C{T)/T is mono- 
tonically decreasing with increasing temperature. The 
tendency towards activated behaviour, signaled by a peak 
at approximately \/ A 2 + e 2 in C(T)/T, is also found at 
all dissipation strengths for sufficiently large (typically 
of order A r ) asymmetries e. In the Kondo language this 
corresponds to the Zeeman splitting of the Kondo res- 
onance due to a local magnetic field. In contrast, the 
dielectric susceptibility, \sb — —d 2 F/de 2 , with F(T) the 
two-level system free energy, shows only a monotonically 
decreasing behaviour with increasing temperature for all 
dissipation strengths < a < 1 in t he sy mmetric case 
(see Fig. ||and for further details Sec. IV C ). 
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FIG. 2. The dielectric susceptibility, \sb, on a linear tem- 
perature scale for the symmetric case at weak and strong dis- 
sipations. 

A peak in Xsb (T) at finite temperature only appears for 
a sufficiently large level asymmetry. The low temperature 
behaviour corresponds to that of a renormalized Fermi 
liquid at all a < 1 with the renormalizations increasing 
with a. Beyond these overall features we also discuss 
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the detailed form of the universal scaling functions of the 
dissipative two-state system for all a < 1 and e. 

The only previous detailed studies of the thermody- 
namics of the Ohmic two-state system which we are aware 
of are, (a), the numerical renormalization group studytll 
and, (b), the work of Gorlich and .Weisala. The latter au- 
thors used a path integral methods to calculate the parti- 
tion function of the dissipative two-state system for both 
Ohmic and non-Ohmic dissipation. Their results are re- 
stricted to weak level asymmetries e <C A r and no results 
are presented for the finite temperature dielectric suscep- 
tibility. For the Ohmic case and a«l they recover the 
linear T behaviour of the specific heat at low temperature 
k B T <C A r and the correct high temperature behaviour 
at all a, however reliable results for strong dissipation 
1/2 < a < 1 could not be obtained within their perturba- 
tive approach. The numerical renormalization group cal- 
culations irJHI are non-perturbative and gave the specific 
heat accurately for all temperatures at both weak and 
strong dissipations. A drawback of—this method, how- 
ever, is the logarithmic discretizationEjtj of the fermionic 
environment. This limits the ability of the method to 
resolve finite temperature features, such as the peak in 
C(T)/T at k B T ~ A r for a < 1. The^lculation of the 
dielectric susceptibility by this methodO is also problem- 
atical at sufficiently low temperatures fcgT < A r i3. As 
explained irJLZI, accurate results for Xsb at T — required 
an analysis of the strong-coupling fixed point Hamil- 
tonian together with the leading irrelevant deviations. 
Thus, this method gave accurate results for Xsb at T = 
(from the fixed point analysis) and for k B T > A r , but it 
was not possible within this method to determine equally 
accurately the bahaviour of \sb for < k B T < A r . As we 
shall see the Bethe Ansatz method we use in this paper 
overcomes all the above difficulties. 

The paper is organized as follows: in Sec. [n] we intro- 
duce the model of the dissipative two-state system and 
outline its equivalence to the anisotropic Kondo model for 
the case of Ohmic dissipation. Some implications of the 
Anderson- Yuval scaling picture of the anisotropic Kondo 
model for the Ohmic two-state system are briefly dis- 
cussed. This gives a qualitative understanding of the 
physics of the latter in both the tunneling and local- 
ized regimes in terms of the fixed points, their stabil- 
ity and their associated low energy scales. Finally we 
show the connection between the scaling picture and the 
renormalization group flow obtained from the exact solu- 
tion of the AKM via the Bethe Ansatz. The correspon- 
dence between the models via bosonization, described 
in Appendix [a|, is then used to translate the thermody- 
namic Bethe Ansatz equations for the anisotropic Kondo 
model, derived by Tsvelik and Wiegman, into the lan- 
guage of the Ohmic two-state system in Sec. [ill] for both 
strong dissipation, a > 1/2, (or weak anisotropy in the 
Kondo model) and weak dissipation, a < 1/2 (or large 
anisotropy in the Kondo model) and at any level asym- 
metry e (or local magnetic field in the Kondo model). 
Analytic results are then presented for the specific heat 



and dielectric susceptibility of the two-state system at 
high and low temperatures and arbitrary dissipation in 
Sec. Ill B 2 and at all tem peratures at the Toulouse point 



(a = 1/2) in Sec. |III B § The Wilson ra tio for the Ohmic 
two-state system is discussed in Sec. IIIB3. Sec. |^ 
gives the numerical solution of the thermodynamic Bethe 
Ansatz equations at all temperatures for both weak and 
strong dissipation and for both symmetric and asymmet- 
ric two-level systems. Our conclusions are summarized 
in Sec. |v|. Appendix |b] contains some details on the 
Bethe Ansatz solution of the AKM and the corresponding 
thermodynamic Bethe Ansatz (TBA) equations which we 
solved in this paper. Appendix |C| gives details of the nu- 
merical procedure used to solve the TBA equations and 
Appendix ^ contains the universal TBA equations for 
weak dissipation (large anisotropies in the AKM), with 
some corrections made to those found originally in Ref. [l3| 



II. MODELS 



A. The dissipative two-state system 



by 



The model of the dissipative two-state system is given 



H S b = ~TiAa x + ^-ea z + '^uj l (a\a i + ]- 



C, 



:(ai+a\). 



(1) 



Here Oi,i = x,y, z are Pauli spin matrices, the two states 
of the system correspond to a z =| and a z —\.. A is the 
bare tunneling matrix element and e is a bias. The en- 
vironment is represented by an infinite set of harmonic 
oscillators (labeled by the index i) with masses m, and 
frequency spectrum Ui coupling linearly to the coordi- 
nate Q — \q a a z of the two-level system via a term char- 
acterized by the couplings Cj. The environment spectral 
function is given in terms of these couplings, oscillator 

c 2 

masses and frequencies by J(uS) = \ * )5(a) — u>i). 

In the case of an Ohmic heat bath, of interest to us 
here, we have J{oS) — 2irauj, for u> << w c , where to c 
is a high energy cut-off and a is a dimensionless pa- 
rameter characterizing the strength of the dissipation. 
The Ohmic two-state model (also called the Ohmic spin- 
boson model) has beea-jntensively studied (for reviews 
we refer the reader tourEl). The model has a low energy 
scale, A r < A for A << uj c , which depends on the 
dissipation strength a, and which may be interpreted 
as a renormalized tunneling amplitude. For a = the 
two-level system is decoupled from the environment and 
A,. = A, whereas with increasing coupling to the en- 
vironment, this energy scale is strongly renormalized: 
A r /uj c ~ (A/tJe) 1 ^ 1- "). Another scale, the frequency of 
tunneling oscillations, Sl(a, A r ) = Q(a)T(a, A r ), is rele- 
vant for time dependent quantities. Here T(a,A r ) ~ A r 
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is the decay rate and Q{a) = cot(^-j-^j) is the qual- 
ity factor of the oscillations^. The latter vanishes at 
the Toulouse point a = i, where the tunneling oscil- 
lations vanish (the "coherence-decoherence" crossover). 
For < a < 1/2 the dynamics coxpesconds to damped 
oscillations of frequency A r )aot3. This is some- 
times called the "coherent"H regime. The system ex- 
hibits phase coherence throught this regime, albeit with 
damped oscillatory contributions to real time dynamical 
quantities. A smooth crossover to "incoherent" behaviour 
occurs at a = 1/2. The tunneling amplitude remains fi- 
nite in the "incoherent" regime 1/2 < a < 1, but there is 
no phase coherence in time dependent dynamical quan- 
tities: f2(a, A r ) = for a > 1/2. Another physically 
relevant value of the dissipation strength is a — 1/3, 
where an inelastic peak, present in the neutron scatter- 
ing cross-section for a < l/Sf-vaaishes and gives rise to a 
quasiel astic pe ak for a > 1 /3c5t3 (see also the discussion 
in Sec. I V B 1 ) . Finally, for sufficiently strong dissipation 
a — ► 1 — , the renormalized tunneling amplitude vanishes 
giving rise to the phenomenon of "localization" or "self- 
trapping" for a > a c ~ 1 (a c depends also on the precise 
value of A). In this paper we will be interested only in 
the thermodynamics of the dissipative two-state system. 
For such quantities the exact solution shows that, in the 
tunneling regime (0 < a < l)c9, the only relevant scale 
is A r . 



B. Equivalence to the anisotropic Kondo model 

The equivalence of the Ohmic two-state system to the 
anisotropic Kondo model (AKM) has been shown at the 
Hamiltonian level via bosonizationu as outlined in Ap- 
pendix [A| This equivalence was believed to be valid in 
the region a > 1/2, which corresponds (see below for 
the precise statement of the equivalence) to the region in 
the parameter space of the AKM between weak-coupling 
(pJiipS^ 1) and the Toulouse point {pJ\\ ~ 1)- Recent 
worked shows that the equivalence extends beyond the 
Toulouse point into the region describing weak dissipa- 
tion < a < 1/2 (or large antiferromaenetic Jy in the 
1, see alsco). The AKM is given bjS 



AKM, 



J\ 



kk' 



~77 X!( c fcT Cfe 'T ~ c ki c k'l)S z + 9HBhS z 



(2) 



kk' 



The first term represents non-interacting conduction elec- 
trons and the second and third terms represent an ex- 
change interaction between a localized spin 1/2 and the 
conduction electrons with strength Jj_, Jy. The last term 
in Eq. (Q) is a local magnetic field, h, coupling only to 
the impurity spin. The correspondence between H and 
Hsb, outlined in Appendix |A|, requires that 



e = g^sh 
A 

— = pJ± 



7T 



(3) 
(4) 

(5) 



where tan 5 = — 7r ^ J|1 , 5 is the phase shift for scatter- 
ing of electrons from a potential Jy/4 and p = 1/2D is 
the conduction electron density of states per spituat the 
Fermi level for a flat band of width 2D = uJ&El. We 
note that weak dissipation (a — > 0) in the Ohmic two- 
state model corresponds to large antiferromagnetic cou- 
pling (Jii — > oo) in the Kondo model whereas dissipation 
strength a > 1 in the Ohmic two-state model corresponds 
to ferromagnetic coupling Jn < in the Kondo model. 



C. Renormalization group flow 

The renormalization group flow of the Ohmic two-state 
system can be obtained by making use of the above Wuiv- 
alence and the Anderson- Yuval scaling equationsES for 
the AKM. These equations hold to lowest order in pj± 
but for all pJ\\: — oo < pj\\ < +oo. They are there- 
fore valid for all < a < 4 provided A/uj c <C 1. The 
Anderson- Yuval scaling equations extend the validity of 
the well known Poor Man's scaling equations to the whole 
Jy axis and reduce to those when pJu < 1. In terms of 
the dimensionless quantities pJ± and 



7T 7T 



(0) 



where S was defj 
equations readc3 

de 



ed above, the Anderson- Yuval scaling 



dlnD 
dpj± 
d\nD 



2)(pJ ± ) 2 +0(pJ ± ) 
0(pJ±_) 3 



(7) 



By using the correspondence between the models given 
above, and noting that e = 2(1 — a), we obtain the fol- 
lowing scaling equations for the Ohmic two-state system 



da 
din ijj c 
d(A/w c ) 
dlrtLj r 



U) r . 



'(1 



a)(£) 



A , 
0{-f 



(8) 
(9) 



Note that in these equations a and A are running vari- 
ables which are functions of the running cut-off, u> c . The 
equations have to be supplemented by specifying initial 
conditions 

a(uj c = uj ) = a Q 
A(uj c = u> ) = A . 
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where ao, Aq, ojq are now the parameters appearing in the 
bare Hamiltonian (where they appeared as a, A, lu c ). We 
shall use this notation for the remainder of this section. 
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FIG. 3. The scaling trajectories of the Ohmic two-state 
system obtained from the Anderson- Yuval scaling equations 
for the AKM. Only the region < a < 3 is shown. The left 
and right separatrices ata = l,A/w c = define the regions 
labeled A, B and C and the arrows indicate the direction of 
decreasing uj c . 



From (|3|-|9|) there is a line of fixed points at A q /lu = 
for ao > 0. Their stability to a finite Aq/ujo follows from 
(^J), which states that A/u> c is relevant, marginal or ir- 
relevant depending on whether the dissipation strength 
ao is less than, equal to or larger than lEd. Hence, the 
line of fixed points at A /ivo = for ao > 1 are stable 
low energy fixed points, whereas the line of fixed points 
at Ao /u)q = for ao < 1 are unstable high energy fixed 
points. The scaling trajectories can be calculated by di- 
viding the two equations (||-||) an d integrating the re- 
sulting equation from o>o down to lu c : 



A. _(*>,> 

LU C UJ 



= — ((In a — a) — (lnao — ao)) (10) 



They are shown in Fig.(|J). The arrows indicate the di- 
rection of decreasing u) c . When the flow is to strong 
coupling, the scaling trajectories will be quantitatively 
correct only for A/lu c <C 1. The scaling diagram is di- 
vided into three regions by two separatrices meeting at 
a = l,A/w c = 0. 

The regime A (ao > 1) corresponds to the localized 
regime of the Ohmic two-state system (or the ferromag- 
netic sector of the AKM). The dimensionless tunneling 
amplitude A/oj c is irrelevant and the flow is to a line of 
fixed points (a = a* and A/lo c = 0). This case is eas- 
ily analyzed since the scaling equations remain valid as 
lo c — > 0. Since A/lu c decreases as uj c — > 0, it follows from 
(||) that a remains almost unrenormalized: a -» a* ~ «o 
as uj c — > 0. Integrating (||) gives a renormalized tunneling 
amplitude A r = A(cj c ) = Ao(w c /wo) Q;o which vanishes at 
T = at low energies. Quantum mechanical tunneling 



is absent for ao > 1 at T = and for sufficiently small 
A /W 
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FIG. 4. Flow of running coupling constants in the localized 
region A of Fig.|^, for ao = 1.4 and Ao/loo = 0.01. 

At finite temperature the low energy cut-off uj c is re- 
placed by ksT resulting in the well known temperature 
temperature dependent tunneling amplitude A r (T) — 
Ao(ksT I 'uJa) ao in the strong dissipation limit. Fig.|] 
shows the flow of the dimensionless coupling constants 
for a typical case in the localized regime. These were 
obtained by integrating (||-^) using the Runge-Kutta al- 
gorithm for 1st order differential equations. 

The regime B (a < 1) corresponds to the tunneling 
regime of the Ohmic two-state system (or the antifer- 
romagnetic sector of the AKM): in this regime A/oj c is 
relevant and the flow for < ao < 1 is away from the 
line of high energy fixed points at A/uj c = towards the 
strong coupling fixed point at a = and A/lo c = oo. 
This is shown in the numerical solution of (^-^), in 
Fig. ||a. The scaling analysis, of course, breaks down 
when A/uj c — O(l), however, other methods, such as the 
numerical renormalization group and the Bethe Ansatz, 
show that the low energy fixed point is at A /cj c = oo and 
a = 0. In this regime, A(lu c ) tends to a finite renormal- 
ized tunneling amplitude, A r as lo c — > 0. In the AKM 
this low energy scale is the Kondo scale, generalized to 
the anisotropic case. We can estimate the T = renor- 
malized tunneling amplitude as the crossover scale sepa- 
rating weak (A/cj c <C 1) and strong coupling (A/u) c ^S> 1) 
regimes of the model. Define A r = A(oj c ) where Cj c is the 
crossover scale such that A(Cj c )/Cj c = 1. Integrating (||) 
down to this crossover scale 



"o d(A/u> c 



A/lu c 



(1 — a)<ilnu; c 



(11) 



and approximating a by ao over this energy range gives 

A r /u; = (AoM,) 1 ^, (12) 

the correct low energy scale for the Ohmic two-state sys- 
tem, up to prefactors depending on ao- 
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FIG. 5. Flow of running coupling constants for (a) «o =0.1 
and Ao/wo = 0.01 corresponding to the tunneling region B 
in Fig. H, and (b) ao = 1.3 and Ao/a;o = 0.4 corresponding 
to the tunneling region C in Fig.^. In (b), A c is a crossover 
scale at which A becomes of O(l). 

Finally, there is a region C, of strong dissipation, which 
also corresponds to the tunneling regime of the Ohmic 
two-state system. However, the flow to the strong cou- 
pling fixed point is such that, for a > 1, A/uj c initially 
decreases with decreasing oj c , signaling a tendency to lo- 
calization for strong dissipation, but eventually, as a re- 
sult of a strong renormalization of a to below 1, A/u> c 
becomes relevant and then increases to strong coupling. 
The flow is to the same strong coupling fixed point as for 
region B and therefore region C belongs to the tunneling 
regime of the model. Due to the strong renormalization 
of a, Fig. ^Jd, it is difficult to estimate the form of the low 
energy scale A r in this case. 

The renormalization group flow described above con- 
sists of a one parameter family of scaling trajectories la- 
beled by a parameter C = C(a, A) which takes a con- 
stant value along each trajectory. This constant of the 
motion is called a scaling invariant and can be found from 

©: 

C(a,A) = -l-i(-) 2 -(lna-a) (13) 
2 uj c 



Corresponding to this one parameter family of scaling 
trajectories we expect the scaling functions for physical 
quantities to consist of a one parameter family labeled 
by C, with different scaling functions for each scaling 
trajectory. The scaling invariant C is not unique but 
depends on the cut-off scheme of the theory. A scaling 
trajectory may be specified differently (i.e. by a different 
function C) depending on the cut-off scheme . One then 
has to identify the scaling trajectories within the differ- 
ent schemes if one wishes to compare results. The Bcthc 
Ansatz solution of the AKM, which we want to use in 
the next section, is an example where such a different 
scheme is used. We state here how we identify the scal- 
ing trajectories of the Bethe Ansatz solution with those 
we discussed above for the Ohmic two-state system (or 
equivalently the AKM with a finite bandwidth 2D = lo c ) 
and leave the details to Appendix |b[ As discussed in the 
appendix, the Bethe Ansatz solution yields a renormal- 
ization group flow depending on two functions, \x and /, 
of the dimensionless couplings of the AKM. The func- 
tion /j, is the scaling invariant and specifies the scaling 
trajectories and the function / sets the low energy scale 



T K = 2Dexp(-f/(jM/n)). 



(14) 



Comparing the high-temperature behaviour of the two 
models we find the correspondence 



fi/w — 1 — OLQ 



(15) 



We shall show in Sec. Ill B 3 that the static susceptibility, 
Xba, calculated from the Bethe Ansatz solution is equal 
to the dielectric susceptibility, \sb, of the Ohmic two- 
state system and that, at T = 0, these are related to the 
Kondo temperature, as defined above, by 



1 



Xba = Xsb 



2w(l- n/n)T K 2-!ra T K ' 



(16) 



This suggests that we define the-cenormalized tunneling 
amplitude A,, in terms of Tk bjc3 



A r = a T u 



K 



(17) 



so that the local T = dielectric susceptibility, Xsb, of 
the Ohmic two-state system is given in terms of A r by: 



Xsb 



1 



2?rA r 



(18) 



The above relations fix / in terms of /i and Tk (or equiv- 
alently in terms of ao and A r ) and will prove useful in 
translating the Bethe Ansatz results of Sec III into results 
for the Ohmic two-state system. 



III. THERMODYNAMIC BETHE ANSATZ 
EQUATIONS 



G 



A. Thermodynamic Bethe Ansatz equations for 

a= 1/v and a = 1 — 1/v 

The thermodynamic Bethe Ansatz equations for the 
anisotropic Kpado model have been derived by Tsvelik 
and WiegmanHI (for a short overview of the derivation see 
Appendix |b|) . For a general anisotropy they consist of an 
infinite set, n = 1,2, . . ., of coupled integral equations for 
the "excitation energies" e„(A) (defined in Appendix ^). 
The thermodynamics is calculated from the impurity free 
energy which depends explicitly only on e\. As discussed 
ir£3, the infinite set of integral equations for the e n de- 
couple to a finite set for values of the anisotropy corre- 
sponding to rational values of the scaling invariant fi/ir of 
the anisotropic Kondo model, and hence to rational val- 
ues of the dissipation strength a = 1 — n/ir in the Ohmic 
two-state system. In particular for anisotropies given by 
fi/ir = 1/v and /i/tt = 1 — 1/v with v = 3,4, . . ., there 
are only v coupled integral equations for the v quantities 
t\, £2, . . . , e u . These two cases allow us to study both the 
weak dissipation (a — 1/v < 1/2) and the strong dissipa- 
tion (a = 1 — 1/v > 1/2) limits of the Ohmic two-state 
system. Explicitly, the v coupled integral equations in 
Eq.6.2.11 ofB are 

^ = s *[ln(l + exp^i)(l + exp^-) 5 



e v -l{X) Wli e ^-2s 

— - — = x + s * ln(l + exp -j^) , 
-jr- =x -s* ln(l + exp -^-) . 



(19) 



Here xo — s^yr, £ is the bias, g = 1 is a ^-factor in the 
Kondo problem and D(X) is the driving term which is 
given explicitly by 

D(X) = — sign[a — ^]"^r arctan(exp(7rA)) . (20) 

with lu c a high energy cut-off (corresponding to the band 
width cutoff, D, in the Kondo model). The integral op- 
erator s* is defined as ir£3: 



a* /(A) 



+ OC 



2cosh(7r(A- A')) 



-f(X') dX'. (21) 



The £j's (j = 1, satisfy boundary conditions at 

A = ±oo that can be obtained easily from Eq. (p9|). The 
impurity contribution to the free energy can be simply 
expressed as 



+oo 



F = -T I s(A + ^-) ln(l + e £l(A)/T ) dA. (22) 

J — oo M 

Eqs. ( |l9| ) and (p2| ) describe the complete thermodynam- 
ics of the model. Note that the field term xq in Eq. jTg| ) 
represents a global magnetic field in the AKM coupling to 



both the impurity and con duction electrons. However, as 
we discuss in Sec. Ill B 3 , due to electron-electron inter- 
actions introduced to assure integrability , the impurity 
susceptibility of the BA solution, —d 2 F/de 2 , coincides 
with the susceptibility of the Ohmic two-state system 
(i.e. with the impurity susceptibility of the AKM with a 
field coupling only to the impurity). Therefore the dielec- 
tric susceptibility, Xsb{T,s), and specific heat, C(T, e), of 
the Ohmic two-state system can be simply calculated as: 



Xsb(T,e) 
C(T,e) 



-d 2 F/de 2 , 
-Td 2 F/dT 2 



(23) 
(24) 



B. Analytic results 

1. Scaling and universality 

A careful analysis of Eqs. (|n]) and ( j22] ) makes im- 
mediately transparent the meaning of universality. For 
fi < tt/2 (a > 1/2) the impurity free energy can be 
shown to be dominated by contributions from the re- 
gion A <C 0. In this limit the driving term can be ap- 
proximated as D(X) w — ¥j?e KX and the explicit cutoff 
dependence can be transformed out of the equations by 
a simple shift, A — > A + -Mn^j-. In this way one ar- 
rives at the following universal equations for the quanti- 
ties^(A)=e 3 -(A+iln^)/T: 

ifj(X) = s* [ln(l + e^- 1 )(l + e^ +1 ) 
+ £ i ,„_ a ln(l + e-*-')] - S lA e* x 
p lJ - 1 (X)=x + s*ln(l + e v "- 2 ) (25) 
</?„(A) = x - s * ln(l + e Vl - 2 ), , 

F = -T f + s(X + - In — ) ln(l + e Vl ^) dX , (26) 

J-oo * T K 

where the Kondo temperature has been introduced as 



T K = w c exp(-7r///i) 
= A r /a . 



(27) 



For /i > 7r/2 (a < 1/2) the seemingly innocent sign 
change of the driving term alters the structure of the 
solutions completely. The derivation of the universal 
equations becomes .much more complicated in this case, 
but is still possibleEJ. These were not used in obtaining 
the numerical results described below, but they served 
as a useful check on the simpler set of equations jlSb . 
We reproduce them, correcting some minor typos inEa, 
in Appendix |d|, together with a comparison of numer- 
ical results obtained from them and the simpler set of 
equations (|l^). The equations above clearly show that 
the thermodynamic quantities depend only on the ratios 
e/T and T/T K - T/A r . Note, however, that while the 
parameter / only influences Tk, for each /x (a) one ob- 
tains a different set of equations and therefore different 
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thermodynamic behaviour. Thus two models have essen- 
tially the same universal behaviour if and only if their 
parameter fi (a) is the same. From these considerations 
immediately follows that the usual RG scaling trajecto- 
ries correspond to the lines Tk = const and p — const. 
(The latter requirement also follows from the fact that 2p 
turns out to be the anomalous dimension characterizing 
the high temperature behaviour, which should be scale 
invariant.) In the small coupling limit one immediately 
obtains the usual leading logarithmic scaling equations 
by expanding / and p in pj± « I± and pJ z ~ I z : 



dpj z 
din uj c 
dpj± 



in agreement with Eq. (Q) 



= -(pJ±f , 
= -p,Jj_p,J z , 



(28) 
(29) 



2. Asymptotic properties 

The asymptotic behaviour of various physical quanti- 
ties can be determined by analyzing Eq. (|l9|). Rewriting 
Eq. (|n]) in terms of the quantities — l,...,v defined 
by 



0=ln[l + exp(^)] (i = l, 
e,=ln[l + exp(-^)] , 



>f-l), 



(30) 



one can easily show that in the A — > — oo limit the asymp- 
totic solution of the BA equations behaves as 



(31) 



with t = 2/i. On the other hand, for A ^> ln(T/w c ) 
£i vanishes extremely fast and can be approximated by 
0. Substituting these into Eq. ( p6| ) and using p/n — 
1 — a and A r ~ Tk one can immediately extract the 
leading behaviour of the impurity free energy. In the 
high temperature limit one obtains: 



F(T » Ar, e) 

/ A rN 2-2a 



T In 



smh(eg/T) 
sinh(e 3 /2T) 

Xs6 (T» A r ; e = 0) = -- 
1 ,Q 2 



8 2 F 



T v 4 
C(T > A r ; e = 0) 



2B 



de 2 

■ A r ^2-2a^ 
• A r ^ 2-2a 



(32) 

(33) 
(34) 



where the constants A and B depend only on p. From 
these equations it is clear that in order to recover the free 
impurity spin at high temperatures — unlike the choice 
g = {y — l)/u of Ref. 13 — we have to take the bare value 



g = 1 for the electronic g-factor. This special choice will 
influence the Wilson ratio discussed below. 

Similarly, for the low temperature regime we obtain: 



(e«T«A,)«-(- + f-^ 
1 9 2 



X sfc (T« A r ; e = 0) 
C(T< A r ; e = Q)- 



A r 2tt 
7r Ta 

3 a7 ' 



(35) 
(36) 
(37) 



where the numerical constants have been calculated fol- 
lowing the same lines as in Ref. Thus at low 
temperatures, the well-known Fermi liquid behaviour is 
recoveredHEj. 



3. Susceptibility and Wilson ratio for the Ohmic two-state 
system 

As discussed by Wiegmann and Tsvelikjil to ensure in- 
tegrability, an artificial electron-electron interaction has 
to be introduced. While this interaction has no effect in 
the course of the solution of the isotropic model, in the 
anisotropical model it renormalizes the electronic and im- 
purity g-f actors: 



<J 



1 



1 



1 



.9 



V 1 - mA 



(38) 



This can be most easily checked by calculating the host 
susceptibility of the AKM , \host = —d 2 Fh os t/dh 2 , fol- 
lowing the lines of Ref. |l3|. (Here Fh os t is the free energy 
of the electrons, and is given by an expression similar to 
Eq. ( p2] ) with / = 0.) After a tedious calculation one 
obtains the result that: 



Xhost 



Xfree 

1 - p/ir 



Xfre 



(39) 



where Xfree = L/(AirN) denotes the free electron sus- 
ceptibility (L is the length of the system, N the number 
of electrons, and vp — fcs = 2/is = 1). Note that the 
specific heat, Ch os t = C/ree = TLtt/N3, is completely 
unaffected by the electron-electron interaction. 



We now prove the statement made in Sec. Ill A, that 



the impurity contribution to the global susceptibility of 
the AKM, obtained from the Bethe Ansatz calculations, 
is identical to the susceptibility of the Ohmic two-state 
system. We denote the bare g- factors of the impurity and 
conduction electrons in the AKM by gi and g e and the 
corresponding susceptibility of the AKM by x(ffi,3e)E3- 
Now in the Bethe Ansatz solution, starting with bare 
values gi = g e = 1, the renormalizations discussed above 
imply that the BA susceptibility, xba = —d 2 F/de 2 , is 
given by 

Xba = x(9i =9e = l/\fa) = x(9i = 9e = l)/a. (40) 
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The dielectric susceptibility of the Ohmic two-state sys- 
tem, Xsb, measures the response to a local electric field 
and is equal to the impurity susceptibility of the AKM, 
xidi = l ; Se = 0), with the magnetic field coupling only 
to the impurity spin: 



Xsb = x(9i = i,g e = o). 



(41) 



This follows from the equivalence of the two models dis- 
cussed in Appendix [A|. We now make use of Eq.21 of 
Rcf. |3^ (valid for arbitrary T), connecting the impurity 
susceptibilities of the AKM with arbitrary g-factors. This 
states that 



X(9i 



1) 



(1 



z-) 2 x(g t 

7T 



l;g e = 0), (42) 



where the phase shifts have been defined in Appendix |A] 
(note the sign chang e with respect to Ref. |36j). Hence, us- 
ing (HIH!) and Eq. (|A16|) a = (1 + 2|) 2 , we find that the 
BA susceptibility is just the susceptibility of the Ohmic 
two-state system, 

X ba(T,s) = -d 2 F/de 2 = X (m = 9e = l/y/a) = Xsb , 

(43) 

as stated earlier. This can be further checked by calcu- 
lating the high temperature Curie susceptibility and the 
Wilson ratio (from Eqs. (|36j) and (E37|)) for which we find 



R s 



Xsb 



Rbj 



X ba(T » T K ~ A r 

Cfree{T) XBA 



lim 



> 4T 
2/a 



(44) 
(45) 



Xfree C (T) 

in agreementp-with exact results obtained for the spin- 
boson model.123 Furthermore, one can easily check by a 
bosonization procedure along the Toulouse lind£j that 
Eq. ( |48|) is identical with the impurity contribution to 
the free energy in the presence of a local magnetic field 
applied at the impurity. We note that the above result for 
the Wilson ratio of the Ohmic two-state system holds for 
all level asymmetries (local magnetic field in the AKM). 
Finally, in order to prevent confusion, we state the con- 
nection of the above Wilson ratio for the Ohmic two- 
state system to that usually encountered in the Kondo 
model. The former is defined using the susceptibility 
of the Ohmic two-state system Xsb which we showed 
was equal to the susceptibility xba resulting from the 
BA calculation on the AKM (with an electron-electron 
interaction to ensure integrability) . The susceptibility 
used in defining the Wilson ratio for magnetic impurities 
is, however, not xba = Xsb = x{9i = l,9e = 0) but 
x(9i = liSe = !)• Therefore, in terms of this susceptibil- 
ity, and using Eq. (ff2|), the corresponding Wilson ratio, 
Rak is given by 

C free (T) x{9i =9e = l) 



Rak 



lim 



(46) 



Xfree C{T) 

The enhancement over the non-interacting value R = 1, 
indicates that the quasiparticles are strongly interacting 
at low temperaturestil. 



4- The Toulouse point: a = 1/2 

The AKM possesses a so-called Toulouse lindHI^II^ 
Along this line the model can be mapped by a simple 
unitary transformation to a resonant level model without 
interaction and can be solved by refermionization. For 
a dissipative two-state system this line has been shown 
to correspond to the special value a = 1/2 (/i = tt/2, 
v = 2 in the BA solution) separating the coherent and 
incoherent tunneling regimes. 

Along this line the BA solution simplifies enormously 
too: For v = 2 only 'one-strings' with parity v = ± 
are allowed, and as one can immediately check by using 
Eq. flBT] ), the rapidities of the spin-excitations are com- 
pletely decoupled. Therefore, in the first of Eqs. ( p~s| ) only 
the driving term remains and one obtains in the scaling 
limit: 

— ei(A) - -— e + — , 



-;£2(AJ - — e + Y 



(47) 



Substituting these expressions into Eq. ( |22|) one immedi- 
ately arrives at 

F = / dk « ln{ 1 + 2 cosh(. 9£ /T) e - fc / T 



(48) 



which coincides with the resonant level result (note the 
slight, difference with respect to the formula 6.2.15 of 
Ref.tj). It is straightforward to verify that in the limit 
T ^ 



Xsb = 9 2 /(ttT k ) = g 2 /2TTA r 



1/2 



and 



C = ttT/(3T k ) = 7rT/6A r , a =1/2 



giving the expected Wilson ratio, R s b = 2/a = 4, for 
the Ohmic two-state system at a = 1/2, with R s b as 
defined above. The high temperature limits S = In 2 and 
Xba{T ^> Tk ~ A r ) = g 2 /AT are also easily verified. 



IV. NUMERICAL RESULTS AT ALL 
TEMPERATURES 

A. Numerical procedure 

The closure of the infinite set of thermodynamic Bethe 
Ansatz equations to a finite set at rational values of the 
dissipation strength a can be used to obtain highly ac- 
curate results for the thermodynamics. This avoids the 
truncation errors associated with solving these equations 
at other values of a. In particular at a = \/v and 
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a = 1 — \jv we have v equations. In the numerical 
procedure we found it more convenient to set up integral 
equations for new quantities — 1, . . . , v defined by 
eq. © 



& =ln[l+exp(^)] = ln[l + «,•], j = l 
&,=ln[l + exp(-^)] =ln[l + « 1/ ] J 



(49) 

^ are introduced for 



where the functions Kj , j — 1 , 
later convenience. 

The TBA equations then take the form, 

&(A) = ln[l + exp(^-,i£>(A) 

+ S * (£,-1 + 0+1 + <*7>-2&>))], 

&/-i(A) = M 1 + exp(x + s * f„_ 2 )], 
0(A) = ln[l + exp(-x + s * £„_ 2 )]. 

The impurity free energy is given by 

F(T, £ ) = -fc s T / + s(A + ///i)a(A, T, e)dA 



(50) 



(51) 



where s(A) = (2 cosh(7rA))~ 1 and ///z is related to the 
low energy scale, T K , of the AKM by Eq. @ 

The Kondo temperature, Tr-, is related to the renor- 
malized tunneling amplitude, A r , by A r = uTk as dis- 
cussed in Sec. He. The entropy, specific heat and dielec- 
tric susceptibility can be obtained by numerically differ- 
entiating the Free energy: 



S(T,e) 
C(T,e) 
x(T,e) 



dF _ 
df ~ 
8 2 F 



-T 



dT^(X,T,e) 
«( A + f/N gf dX 

d 2 T&(A,T,e) 



dT 2 



d 2 F 
W 2 



4T 



S (A + 1 lii) 



or 2 

d 2 ei(A,T, £ ) 
dx 2 



dX 



dX 



A more accurate procedure is to set up integral equations 
for the derivatives d (T£j)/dT, d^j/dx and d 2 £j/dx 2 ] . 
More precisely, we set up integral equations for a new 
set of functions, Ej,Fj and Gj = d (TQ/dT = & + 
Td^j/dT where Ej and Fj are the first and second field 
derivatives of 



F, 



dx 



dx 2 



dnj/dxo 

1 + Kj 

d 2 Kj/dxl 

1 + Kj 



diij I Oxq 



and the functions k 9 - where defined in (kt£ 



(52) 



(53) 



Each set 



of functions Ej : Fj,Gj,j — 1, . . . , v then obey coupled 
linear inhomogeneous integral equations. The equations 
are, 

Ej = (1 - e-e>) 8 * {Ej-i + E j+ i + S j>v - 2 E V ) 
^_ 1 = (l-e-«"- 1 )(s*^_ 2 + l) ! (54) 
E v = (1 - e -6 ') (s * E v - 2 ~ 1) 



for the Ej , with the inhomogeneity appearing in the last 
two equations. For the Fj we have 

Fj = Qj + (1 - e -fc) s * {Fj_ x + F j+1 + S itV - 2 F v ) 



(55) 



F u -i = Q.-i + (1 - e"^- 1 ) a * F„_ 2 , 
F„ = Q i/ + (l- e -«-)s*F l/ _ 2 , 

with an inhomogeneous term, 

Qj =E 2 e~^(l -e-**)- 1 . 

and for the Gj we find 

Gj = Sj + (1 - e-**) s * (Gj-! + G j+ i + Sj, u - 2 G V ) 
G v -x =5„_i + (l-e-^- 1 )s*G t/ _ 2 , (56) 
G v = S v + (l-e-Z»)s*G v - 2 

with an inhomogeneous term, 



-e ^'lne ^ — (1 — e ^)ln(l — e 



The numerical procedure is to first solve iteratively 
for the = 1, ...,v. Using this solution one then 

iteratively solves for the functions Ej, Fj , Gj, j — 1 , . . . , v 
in turn. Fig.^ shows a graphical representation of these 
integral equations. 



1 

o- 



2 

o 



3 

o 




v-1 

FIG. 6. Graphical representation of the integral equations 
for the t^E^Fj^jJ = l,...,u. 

The entropy and susceptibility are thereby obtained 
without the need to take any numerical derivatives. Only 
one derivative is required to obtain the specific heat from 
the entropy, so we did not set up separate integral equa- 
tions for the second temperature derivative of T£i . Such 
a procedur£_has been used for zero field static suscep- 
tibilities inE3. This approach also overcomes difficulties 
at large fields (level asymmetries) and low temperatures 
found in early .treatments of similar TBA equations for 
Kondo systemsa provided one deals with the exponen- 
tial decrease in the j = 1 (j = v) functions for strong 
(weak) dissipation respectively. This and further details 
of the numerical procedure and its accuracy are given in 
Appendix ^|. 
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1. Choice of parameters 

The TBA equations were solved for weak dissipation 
at a = 1/6,1/5,1/4,1/3, and for strong dissipation at 
a = 2/3,3/4,4/5. The exact closed solution was used 
to obtain results at the Toulouse point a = 1/2. The 
thermodynamics was calculated at temperatures t m = 
ak B T m /A r = 2" 1 / 2 with m = -20, -19, . . . , +19, +20, 
and for level asymmetries e n = ae n /A r = 2" with n = 
—4, —3, . . . , +3, +4 and for the symmetric case e = 0. 



B. Entropy and Specific heat 

1. Symmetric case: e = 

The entropy of the symmetric two-state system is 
shown in Fig.|?]a as a function of temperature for several 
values of the dimensionless dissipation strength, a, rang- 
ing from weak to strong dissipation. The correct value of 
the entropy, S = In 2, is recovered at high temperature 
for all a. Fig.^b shows the universal specific heat curves 
for the dissipative two-state system. As in other strongly 
correlated systems^ we observe a characteristic cross- 
ing point for the specific heat curves at a temperature 
ksT + /A,. = 0.66 + 0.02. At low temperature the specific 
heat is given by C{T) = aj{T/A r ) + 6(a) (T/A r ) 3 + . . . 
with a linear coefficient of specific heat 7 = aj / A r which 
vanishes as a — > 0. From the definition of the low energy 
scale A r in terms of the zero temperature susceptibility 
and the Wilson ratio (to be discussed below) it follows 
that 7 = 7r/3 for all a, a useful check on the numerical 
results. The coefficient b(a) of the T 3 term is negative 
for a > 1/3, a special point in the parameter space of the 
dissipative two-state system. 
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FIG. 7. (a) Entropy, S(T), and (b) specific heats, C(T)/T, 
for the symmetric two-state system (e = 0) for weak (a < 1/2) 
and strong (a > 1/2) dissipation cases. The T 3 coefficient m 
C(T)/T is negative for a > 1/3 and positive for a < 1/3. 

The significance of a = 1/3 is best seen in the context 
of the dynamics of the two-state system, where it corre- 
sponds to the value of the dissipation strength at which 
the frequency, £l(a), of tunneling oscillations (manifested 
in real time correlation functions) becomes equal to the 
decay rate, T(a), of these oscillations, i.e. fi(a = 1/3) = 
T(a = 1/3) (ojhrthe quality factor Q(a) = O(a) /T(a) 
becomes unity )qcJ. For dissipation a < 1/3 we have 
il(a) > r(a) and the well defined oscillatory mode ap- 
pears to be reflected— in the characteristic peak in the 
specific heat C(T) /Tel 
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a = 2/3 
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2. Asymmetric case: e > 
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FIG. 8. (a) Specific heat, C(T)/(aT/A r ), showing the 
development of the peak at fcsT w A r for a < 1/3 and 
e = 0. liniT^o C(T)/(aT/A r ) = 7 = tt/3 = 1.04719755 is 
recovered to 5 decimal places, (b) the logarithmic derivative, 
dlog C(T, e)/dlog T, at Q = 2/3 for the symmetric (e = 0) 
and asymmetric (e/A r = 8) cases. The approach to the ex- 
pected power law C(T) ~ (A r /k B T) s with 5 = 2a - 2 = 2/3 
for a = 2/3 is found at high temperatures for both cases. For 
the asymmetric system the power law arises at higher temper- 
ature corresponding to the higher low energy scale behaving 
as Ve 2 + A? > A r for e > A r . 

The peak in C(T)/T for a < 1/3 is shown in more 
detail in Fig.||a and is reminiscent of the activated be- 
haviour seen in non-interacting two-level systems. Since 
the excitation spectrum of the Ohmic two-state system is 
gapless for all a > 0, there is no exponential suppression 
of C{T)/T at low ksT < A r as with non-interacting 
two-level systems. The linear specific heat persists at 
low temperature and the system remains strongly inter- 
acting down to T — as is also clear from the value 
of the Wilson ratio (see later). At high temperatures 
ksT ^> A r the specific heat vanishes according to a 
power law, C(T) ~ (A r /fcsT) 2_2a with an a depen- 
dent power in accordance with the asymptotic result 
given by Eq.(|34|) and inlia. This is shown in Fig.||b and 
we see that, as for low temperatures, the behaviour at 
high temperatures is again drastically different to the 
behaviour of a non-interacting two-level system which 
shows C(T) ~ (Ao/fcsT) 2 for k B T > A with A the 
bare tunneling matrix element. The limit a — ► 1~ corre- 
sponds to the weak coupling Kondo model. For weak cou- 
pling, pJ±, pJu <C 1, the exchange operators in the AKM 
become marginal (see above discussion of the Anderson- 
Yuval scaling equations) and consequently the high tem- 
perature properties acquire logarithmic corrections lead- 
ing to a slow approach of the thermodynamic quantities 
to their high temperature limitalJ ( this can be seen in 
the entropy in Fig. ^|a and more clearly in the results for 
the susceptibility to be described in the next section). 
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FIG. 9. Entropy, S(T), for the asymmetric two-state 
system for a range of asymmetries at (a) a — 1/5 and (b) 
a = 4/5. 

We now turn to the asymmetric two-state system. 
Fig. [)] shows the temperature dependence of the entropy 
for different level asymmetries at a = 1/5 and a = 4/5. 
The correct high temperature limit S = In 2 is recov- 
ered for all level asymmetries and dissipation strengths. 
As for the symmetric case, we see again that the en- 
tropy approaches its high temperature limit more slowly 
for strong dissipation than for weak dissipation, again a 
result of increasing marginality of the interactions with 
increasing a at the high energy fixed point. 

In Fig. p~0| — [Tl] we show the specific heats for different 
level asymmetries and dissipation strengths. The specific 
heat remains linear at low temperature, C(T, e) ~ 7T, for 
all level asymmetries. The linear coefficient 7 ~ 1/A r is 
reduced with increasing asymmetry e, a consequence of 
the increasing low energy scale with increasing e, A r — > 
V^A 2 , + e 2 . We see that a sufficiently large asymmetry 
eventually leads to a peak in C(T)/T for all dissipation 
strengths, but that for such a peak to form requires a 
sizeable asymmetry for a > 1/3. 
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FIG. 10. Specific heats, C(T)/T, for the asymmetric 
two-state system for a range of asymmetries, and some typical 
cases for weak dissipation, (a) a = 1/5, (b) a — 1/3. 
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FIG. 11. Specific heats, C(T)/T, for the asymmetric 
two-state system for a range of asymmetries, and some typical 
cases for strong dissipation, (a) a = 2/3, (b) a — 4/5. 



It is important to note that even for large asymme- 
tries, the shape of the specific heat curves is still dif- 
ferent to those of a non-interacting two-level system: at 
low temperature the specific heat remains linear rather 
than exponential and at high temperature the asymp- 
totic behaviour of C(T) is not the non-interacting 1/T 2 
but instead behaves as \/T 2 ~ 2a as shown analytically for 
e = and numerically for s > in Fig. ^| Only for a < 1 
do we expect the specific heat to be reasonably described 
by the non-interacting result, and then only outside the 
Fermi liquid regime fc^T > A r . 



C. Dielectric Susceptibility 

1. Symmetric case: e — 

Fig. [l2] shows that the dielectric susceptibility of the 
dissipative two-state system remains finite down to T = 
for all dissipation strengths a < 1. This is also shown 
in Fig. |l3] together with the Fermi liquid T 2 corrections 
at k B T < A r given by Eq.©, X sb{T) = Xrf(Q)(l " 
c{a)(k B T/A r ) 2 ). By our definition X sb{T = 0) = 
l/27rA r , we have that A r x s f,(T = 0) = l/2vr = 0.1591549 
which is reproduced by our numerical solution to 5 dec- 
imal places in all cases (Fig. |l^). In contrast to the 
specific heat, C(T)/T, the susceptibility is a monotoni- 
cally decreasing function of temperature for all dissipa- 
tion strengths. There is no signature of the onset of ac- 
tivated behaviour in the susceptibility as there was for 
a < 1/3 in C(T)/T. As we shall see below, a finite 
temperature peak in \sb only arises when there is a fi- 
nite level asymmetry. The dielectric susceptibility looks, 
superficially, like that for a non-interacting system, how- 
ever the universal scaling curves depend sensitively on a, 
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as can be seen in Fig. [l^, so that this resemblance is 
misleading. 



0.20 




FIG. 12. Dielectric susceptibility, X*>b(T), for the symmet- 
ric two-state system (e = 0) for weak (a < 1/2) and strong 
(a > 1/2) dissipation cases. The susceptibility is finite at 
T = with A r (a)x(T = 0) = 1/2-7T in all cases as seen in (a) 
and attains its free-spin value of 1/4T at high temperatures 
fcflT S> A r , as seen in (b). 

The strong renormalization of the tunneling amplitude 



A r /uj c ~ (A /w t 



\l/(l-a) 



as a — > 1 gives rise to strongly 



renormalized dielectric susceptibilities at low tempera- 
tures and strong dissipation (xsb(T = 0) = l/27rA r ). 
The approach of the susceptibility to its free spin value 
of 1/4 at k B T > A r at high temperatures (Fig.pjb 
) is governed by power laws with exponents which are 
functions of the dissipation strength as given by Eq.(33), 
k B Txs b {T) « 1/4- 2B{A r /k B T) 2 - 2a , and verified in 
Fig. [ill The approach to the free spin value becomes 
slower as a — > 1~ and eventually logarithmic corrections 
to the susceptibility set in (see Fig.[12|b ). 



e 

X 



0.1592 



0.1590 



0.1588 



0.1586 




0.000 



0.001 

(k B T/A r ) 2 



0.002 



FIG. 13. The T 2 Fermi liquid corrections to the dielectric 
susceptibility at low temperatures for the symmetric case and 
q = 1/6, 1/5,... ,3/4, 4/5. 




FIG. 14. High temperature corrections to the dielec- 
tric susceptibility for a = 2/3 for e — and e = 8. 
lim^oo log 10 [(fc fl T) 2Q - 2 (l/4 - k B T X sb(T))] = log 10 (2B) 



2. Asymmetric case: e > 



The dielectric susceptibility in the presence of a level 
asymmetry is shown in Fig. |i~5|-|T6|. For all dissipation 
strengths we see that a sizeable asymmetry of the order 
of A r is required to give a finite temperature peak in Xsb- 
The approach of the susceptibility to its high tempera- 
ture limit of 1 /AT is governed by the same power laws as 
those found for the symmetric case and verified in Fig.|l4|. 
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FIG. 15. Dielectric susceptibility, Xsb{T), for a range of 
asymmetries and some typical cases of weak dissipation, (a) 
a = 1/5, (b) a = 1/3. 



0.20 



0.15 



0.10 



0.05 - 



0.00 



0.20 



0.15 - 




a = 4/5 



0.10 



0.05 



0.00 




10 10 



10 10 10 10 



kJ/A 



FIG. 16. Dielectric susceptibility, Xab(T), for a range of 
asymmetries and some typical cases of strong dissipation, (a) 
a = 2/3, (b) a = 4/5. 



3. Wilson ratio 



V. SUMMARY AND DISCUSSION 



The Wilson ratio for the Ohmic two-state model, R s b, 
was defined earlier together with the usual Wilson ratio, 
Rakm, for the AKM. These take the values 2/ a and 2 
respectively and are valid for both the symmetric and 
asymmetric cases (i.e. in both zero and finite fields for 
the Kondo model). The Wilson ratio served as a useful 
check on our numerical solution, which recovered it with 
an accuracy of not less than 4 decimal places for all a 
and e. 



In the present paper we studied the thermodynamics 
of a TSS with Ohmic dissipation by exploiting a mapping 
between the DTSS and the anisotropic Kondo model, and 
solving the Bethe Ansatz equations derived by Wiegman 
and Tsvelik for the latter. Treating the BA equations in 
a careful way we were able to calculate essentially exactly 
the specific heat and the susceptibility of the DTSS for 
all temperatures and level asymmetries in the delocalized 
phase of the DTSS, characterized by dissipation strengths 
< a < 1. The Bethe Ansatz solution makes the univer- 
sal properties of the DTSS clear: thermodynamic quan- 
tities are universal functions of two variables, fcsT/A r 
and e/ksT for all a < 1. In the limit a — > 1~ these 
functions reduce to those of the usual isotropic S = 1/2 
Kondo model, which is seen .as a special point in the 
parameter space of the DTSSU. The well known loga- 
rithmic corrections to physical quantities for a — ► 1~ at 
high temperatures fcgT 3> A r give way to power law cor- 
rections away from a = 1. We determined these power 
laws both analytically and numerically for all a < 1 at 
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finite e. In the context of the RG the change from log- 
arithmic to power law corrections to physical quantities 
at high temperature indicates that the tunneling term 
in the Hamiltonian changes from being marginally rel- 
evant (about the high energy fixed p oint) at a = 1 to 
relevant at a < 1 as discussed in Sec. 1IC. At low tem- 



perature ksT <C A r the thermodynamics of the DTSS 
is that of a renormalized local Fermi liquid with an en- 
hanced linear specific heat C(T) ~ akBT/A r , and an 
enhanced, but finite, dielectric susceptibility at T = 0, 
Xsb{T = 0) = l/2nA r . The renormalizations increase 
dramatically as a — > 1~ due to the strong renormaliza- 
tion of the low energy scale A r /uj c ~ (A/uj c ) 1 ^ 1 ^ a ' > . 

We have shown that the characteristic thermodynamic 
properties of the DTSS change smoothly as one increases 
the dissipation strength from a <C 1/2 to 1 > a > 1/2, 
corresponding to weak and strong dissipations, respec- 
tively. In the former case, where the DTSS displays co- 
herent oscillations between its two positions^, we find 
the expected tendency towards activated behaviour in 
the specific heat. A clear signal of this behaviour is 
the appearance of a peak in C(T)/T at the renormal- 
ized tunneling amplitude A r . Such a peak is absent for 
dissipations a > 1/3 in the symmetric case. A finite level 
asymmetry accentuates the tendency towards activated 
behaviour and always gives rise to a finite temperature 
peak in C(T)/T at temperature k B T = A r ~ \/A% + e 2 
provided e > A r . For strong dissipation a > 1/2, the 
specific heat is qualitatively similar to that of isotropic 
5 = 1/2 Kondo systems and shows a monotonically de- 
creasing C(T)/T with increasing temperature. The di- 
electric susceptibility of the DTSS was calculated at all 
dissipation strengths, level asymmetries and tempera- 
tures for the first time. For the symmetric case we found 
that this quantity decreases monotonically with increas- 
ing temperature for all dissipation strengths and develops 
a finite temperature peak only for sufficiently large level 
asymmetries (of the order of A r ). 

The Ohmic two-state system is a generic model capable 
of describing a large number of different physical systems. 
Previous theoretical and experimental work on such sys- 
tems has, however, largely focussed on the dynamic prop- 
erties. Given the detailed understanding which we now 
have of the thermodynamics, it may be worthwhile to 
consider also thermodynamic measurements on Ohmic 
two-state systems. We therefore briefly mention below 
some possibilities where our results could be directly 
tested. 

One of the possible physical realizations is provided by 
two-level systems in metals. Amorphous metals are not 
the best candidates, however, since in these materials a 
broad distribution of DTSS's with very different physical 
parameters occur. Therefore, to calculate their contribu- 
tion to the thermodynamic properties one should average 
over them. This averaging is not impossible, but would 
require solving the BA equations for arbitrary values of 
a, something whiGh-,could be implemented by using fur- 
ther results of Refill In addition one would require the 



form of the TSS-distribution. Better candidates are H 
tunneling in Nbo and metallic materials with tunnel- 
ing centers focmed by substitutional impurities, such as 
Pbi~ x Ge x Tet3. Good quality single crystals can be pro- 
duced from the latter alloy. The Germanium ions form 
identical eight-state systems, which couple to the con- 
duction electrons. Applying external stress on the sam- 
ple one can reduce the degeneracy of the lowest lying 
states to two. Since in this case the DTSS's have approxi- 
mately identical parameters for weak concentrations their 
individual thermodynamic properties can be observed by 
measuring the sample's macroscopic properties. 

Another possible candidate for thermodynamic mea- 
surements is provided by a SQUID. In this case the en- 
ergy difference between the two flux states of the SQUID 
can be easily tuned by an external magnetic field. Mea- 
suring the average flux ($) as a function of the asymme- 
try energy one can readily determine the susceptibility 
X ~ d($>)/dB that can directly be compared to our cal- 
culations. 

Finally, the results we have obtained for the thermody- 
namics of the anisotropic Kondo model may have some 
relevance to the problem of an isotropic S = 1/2 Konda 
impurity in a Luttinger liquid. Schiller and IngersentEj 
have shown that a S = 1/2 Kondo impurity in a modi- 
fied Luttinger liquid consisting only of left-moving spin- 
down electrons and right-moving spin-up electrons can 
be mapped exactly onto the AKM. It is clear that in a 
more realistic description of the Luttinger liquid that this 
exact mapping will only be approximate, nevertheless we 
expect that some of the general trends in the thermo- 
dynamic properties of a S = 1/2 Kondo impurity in an 
interacting system should be captured by such a mapping 
onto the AKM or equivalently onto the Ohmic two-state 
system. The strength, U, of the Coulomb interaction 
can then be relatedE 2 ] to Jii in the AKM and hence to the 
dissipation strength a in the Ohmic two-state system. 
The non-interacting system with U = corresponds to 
a = 1 and U — > oo corresponds to a = 0. Within such 
a picture we therefore expect, from the results of this 
paper, that the thermodynamic scaling functions of an 
isotropic S — 1/2 Kondo impurity in a Luttinger liquid 
will change continuously with increasing Coulomb inter- 
action U (or decreasing dissipation strength a). At low 
temperature the thermodynamics of such systems should 
be similar to, that of a local Fermi liquid, as recently 
discussed ir£3. For sufficiently large U (i.e. sufficiently 
small a), C(T) /T may develop a finite temperature peak 
at ksT ~ A r ~ Tk in analogy to our findings for the 
Ohmic two-state system. Such a peak could then be 
taken as a signature of the Kondo effect in a strongly in- 
teracting system. The unusual heavy fermion behaviour 
of Nd2-pCe x Cu04, e.g. the non-monotonic behaviour of 
C(T)/TO, may be consistent with such an interpreta- 
tion. Strong interactions have also been invoked to ex- 
plain this behaviour in Ref.Lj. Other, more conventional 
explanations, such as lattice coherence effects, cannot be 
ruled out however. 
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The approach we have developed in this paper can be 
extended to more complicated models of two-level sys- 
tems coupled Ohmically to an environment. Thfi-jsffect 
of indirect or electron-assisted tunneling processesEKll on 
the thermodynamics of two-level systems in metals will 
be studied in a future publication. 
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APPENDIX A: EQUIVALENCE OF THE AKM 
TO THE OHMIC TWO-STATE SYSTEM 

1. Bosonization of free fermions 



(A4) 



Consider first the free fcrmion Hamiltonian, 
H Q = J2 £kC U Ck ^- 

k ji 



(Al) 



We take a linear dispersion relation for the conduction 
electrons, Sk = v F k, and measure k relative to the Fermi 
wavenumber kp. We introduce fermion fields ipp,(x) via 
a Fourier series, 



^(x)=L- 1 / 2 J2e~ ik3: c ktt 



(A2) 



with wavenumbers k = 2Tm/L,n — 0, ±1, ... for peri- 
odic boundary conditions appropriate to a finite system 
of length L. The density of states per spin direction is 
given by p — \j2-KVp. The kinetic energy can then be 
written as, 

r+L/2 

H = v F Vk[ c fe(i = iv F / ipj l (x)d x tp ft (x) dx 
k,n J - L ' 2 

(A3) 

The fields Vv( x ) are expressed in terms of hermitian 
bosonic fields tp(x) in the standard way, 



where a is a lattice spacing, required for obtaining con- 
vergent momentum sums. The are the Klein factors 
required to ladder between states with different fermion 
number, and to ensure the carrect anticommutation re- 
lations for the fermion fieldsEjtj, and 



with the bosonic fields <f>, <jr given by 



4>Ux) = (^(a;))t 



J2< 



/2_igx t p - a <l/ 2 



(A5) 



(A6) 



q>0 



The a q , a q , defined for q > 0, satisfy boson commutation 
relations with n q = (qL/2Tr) 1 / 2 , and are given by, 



,t = 



t _ 



-i/2 r 1 t 



(A7) 



It is convenient to introduce spin and charge density op- 
erators in place of a q -\ , a q [ as follows: 

1 / 

a qC = -j= (a gT +a qi ) 

a qS = -j= Kt - a qi) ■ 
The corresponding hermitian fields are, 

with commutation relations, 

[(pc(x), <fis(x')} = 0, 
[<Pc,s( x )> dx"Pc,s(rf)] = 2ni5{x - x 1 ) 



(A8) 
(A9) 



In terms of these we have for the fields and densities, 



'2ira 



(A10) 



Pc,si x ) = -^=^ cs {x)^ c ,s(x) = ~^75 dxlf c,s(x) (All) 



2. Transformation to the Ohmic spin-boson model 

We start with the anisotropic Kondo model 

H akm = H +H ± + (A12) 
with Ho as above and, 
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H± = -y ( c U c k'l S + 4i c fc'T S+ ) 
k,k' 

H W = ~2 E ^T Cfe 'T " c ki c k'l)S z 

k,k' 

In terms of <pc,S, we can write, 
H = v F ^2 1 (4t a «T + a li a il ) 



/+L/2 j 
— : (d x (p c (x)f + (d x <p s {x)f : 
-L/2 27r 

= ^5,(^(0)V T (0)- Vl (0)^(0)) 

= ^(V}(0)V i (0)5- + Vl(0)^(0)5 t ) 

= A C e i^.(o) f f t 5 - + e -i^.(o) f f 4 5 +) 
47ra V 1 i J 

We note that tpc (which commutes with <pg) does not 
couple to the impurity and only gives a contribution to 
the kinetic energy. 



H 



Hj 



b. Identification of the parameters of the Kondo 
Hamiltonian with those of the Spin-boson Hamiltonian 

We notice that a q c does not couple to the impurity 
so we write a q = ia q s and omit the charge density op- 
erators to obtain the Hamiltonian for the spin density 
excitations, 



Hsi 



-7Z-^ + v F ^2qa\a q 



Ana 



J|| 
4tt 



VF 



<?>o 



This is precisely the spin-boson model, 



50 



J q "-q^q 



E 



Co 



q>0 



■at) e^l\ 



[a q + a' q ) 



with oj q = v F q and a spectral function for the harmonic 
oscillators, 



2 '■ — ' m n uj, 
q 



q^q 



2ira 

— 2~~ ' 

% 



a. Canonical transformation on the bosonized Model 

We show that the unitary transformation U = 
exp(i\/2S z ips(0)) applied to Hakm gives the spin-boson 
Hamiltonian, Hsb, for Ohmic dissipation, i.e. that 
V HakmU^ — Hsb- We use the Baker-Hausdorff for- 
mula e~ B Ae B = A + [A, B] with [A, B] a c-number and 
the commutation relations for (pc,tpSi to obtain, 



UH U^ =H - V2v f S, 



dips 



dx 



x=0 



H - V2vfS z 

T 



q 

UH\\U^ =H\\+ constant, 



Ana 1 1 



On using the identities, 



-iV2tp s (0) 



s~ 



US+U^ = e ^s(o)5+ ; 

and the representation, \a + = F^F^ S + , o~ = (er+)t 
— S z , of the Pauli spin operators, the term 



and i<7. 



UH±W becomes, 



UH ± U* = J ± -Lct x . 

4ira 



(A13) 



provided one chooses, 



C a 



m a 



-2 f 2nv F \ 



1/2 



UJ a e 



with a cutoff, 



a 



One can identify the parameters, 



A 
~2 



J± 



and 



4irv F 



which, together with the density of states (per spin) of 
the conduction electrons, p = 1/2t:v f , and the above 
definition of lo c , result in the following identification of 
the dimensionless couplings of the two models, 



A 

— = -pJ± 



(A14) 
(A15) 



The sign of J± (or A) plays no role and we may choose 
A/ui c = +pj ± > 0. 

We remark here that the precise form of Eq. ( Al5| ) de- 
pends on the specific regularization scheme used. Within 
the framework of Abelian bosonization the coupling J» is 
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directly proportional to the phase shifts, pj\\ — An/ 5 
For a finite-band model with cut-off lu c , the expres- 
sion for a in terms of pJ\\ the p hase shift is given by 
5 = — arctan(7rpJ||/4), and Eq. ( A15 ) takes the more 
general form: 



(1 



(A16) 



APPENDIX B: DERIVATION OF THE 
THERMODYNAMIC AL BETHE ANSATZ 
EQUATIONS 

1. Bethe Ansatz equations for finite system sizes 

To construct the Bethe Ansatz (BA) solution of the 
model one first rewrites the AKM Hamiltonian Eq. (||) 
in coordinate space as: 

N r 



H = ^{-idj + ^6( Xj )[l z a?S z + + h.c)] } , 

(Bl) 



where Xj and dj denote the coordinate and spin of con- 
duction electron j, (J = 1, ..,N e ), and S = do/2 is the 
impurity spin. The explicit relationship between the di- 
mensionless couplings I Zt ± and the J's in Eq. (|J) is not 
straightforward and depends on the specific regulariza- 
tion and cutoff scheme used in the original model and the 
BA approach. For small couplings I z « J z and I± ~ J±. 

As a next step one constructs the impurity-conduction 
electron scattering matrix from Eq. (Bl). Introducing 
the impurity and the conduction electron 'rapidities', 
Ao = — 1 and Xj — 0, the s-matrix between the j'th 
electron and the local impurity spin can be written as: 



(B2) 



where — apart from some unimportant phase factors - 
Soj can be expressed as 

S oj (X) = a(A)P TT + b(X)P n + \c(\)(a+ <g> crj + h.c.) 



a(X) sinh(i/i + A/) 
c(A) sinh(i/i) 
a(A) sinh(i/i + A/) 
b(X) = sinh(A/) 



(B3) 



with Pjj and P^ being the projection operators for par- 
allel and opposite spins. 

The parameters \x and / are connected to the bare 
couplings via: 



cos(^i) 



cos(7_l/2) ! 



cth 2 (f) = 



3 in 2 (/||/2) 



3in((7||+7 ± )/2)sin((7||-7 ± )/2) 



(B4) 
(B5) 



The electron-electron s-matrix is not determined by 
Eq. (Bl) but is fixed by the requirement of integrability, 
and is simply given by Sij(Xi — Xj). With this choice the 
impurity-electron and electron-electron s-m atrices belong 
to the same family of s-matrices, Eq. (B3), and satisfy 
the following Yang-Baxter relations: 

S l3 (a)S, k {a + f3)S jk {0) = S jk (j3)S ik (a + fflSy (a) 

(B6) 

insuring the integrability of the modelEl'S. The Bethe 
Ansatz (BA) equations are derived using standard alge- 
braic methods, based on Eq. (B6) by creatiHg-^pin waves' 
from the completely spin-polarized statetL3'E2l. The spin 
waves are characterized by rapidities A/j (/3 = 1,..,M) 
and the obtained state is an eigenstate of H if these ra- 
pidities satisfy the following BA equations: 



M 

n 



sinh(/i(A Q — Xp + i)) _ sinh(^(A Q + f/fi + i/2)) 
-J- x sinh(^(A Q - A/3 - i)) sinh(^(A Q + f/fi - i/2)) 

nJV e 



sinh(/i(A Q + i/2)) 



sinh(/x(Ao, — i/2)) 



(B7) 



Thus an eigenstate of the Hamiltonian is characterized by 
the momenta kj (j = 1..N) of the conduction electrons 
and the rapidities {A^}, corresponding to the charge and 
spin degrees of freedom, respectively. The energy of a 
state is determined through the periodic boundary con- 
ditions: 



M 

n 



sinh(^(A /3 +i/2)) 
^ sinh(^(A /3 - i/2)) 1 



N 



(B8) 
(B9) 



E = Y,kj-hgS^ . 

Since the X^s (/3 = 1, .., M) are spin —1 excitationsB the 
total spin is simply S z tot) = = (N e + 1 — 2M)/2. 



2. Thermodynamic BA equations 

To derive the thermodynamic BA equations (TBA) one 
takes the limit L,N — > oo, N/L = D. In this limit the 
rapidities A Q in Eq. (B7) can be shown to organize into 



{r.v) 



(/3 = 1, .., M/ r<v -\), with integer 



r-stringgi3, A Q — > A 
ranks r and parity v = ±: 

X (r,v) ^ A (r,«) ^ = A („,) + [r±l _ ? ] +»*_(!_ v ) . 

(BIO) 

The allowed (r, v) values must satisfy the stability con- 
dition 

usin(M)sin(/z(r-g)) >0 , q = 1, .., [r/2] , (Bll) 
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and can be determined by expressing /u/tt as an infinite 
(finite) fraction. For rational /z/tt's oaiy a finite number 
of possible ranks and parities survivecil. Specifically, for 
H = -k/v only the pairs 

n = (r, v) = (1, +), (2, +), ..,(1/- f , +) and (1, -) (Bf 2) 

are allowed, which we label by n — 1, .., v in what follows. 

After classifying the possible strings one introduces the 
density of rapidities of the n = (r, u)-strings, g n (X), and 
the density of 'holes', g n (X). These are connected by a 
linear integral equation derived from Eq. (B7) as 



g n (X) = a n (X) -J2[ d\'K nn ,{\ - \') Bn ,{\') , (B13) 

n' 

where the functions a n and the Kernel K nn i depend only 
on fi. The spin part of the energy, Eq. (BE), and the 
entropy can be expressed as: 

E (spin) = J2 J dXf n (X)g n (X) , 



(B14) 



S (s P in) =k B J2 d\[ Qn \n(l + ^) + g n bx(l + ^)] 

„ J Qn Qn 



Then Eqs. ( |19| ) and (g2|) can be derived simply by mini- 
mizing the free energy F^™) = - TS {spm ^ with 
respect to <5g„(A) using the constraint Eq. ( |B13| ) and in- 
troducing the 'excitation energies' e„(A)/T = ln(g n /g n ). 



APPENDIX C: NUMERICAL SOLUTION THE 
THE THERMODYNAMIC BETHE ANSATZ 
EQUATIONS 

In this appendix we describe in detail a new numeri- 
cal procedure for solving thermodynamic Bethe Ansatz 
equations of the form ( |5(i| ) for the quantities = 
Such equations arise also in other contexts, 
for example- in the context of the anisotropic Heisen- 
berg modeliJ. The method we describe has the advan- 
tage that the thermodynamics is calculated without the 
need to take numerical derivatives of the free energy 
(see Sec.IVa). An implementation of this for the spe- 
cial, case of zero field static susceptibilities has been given 
ir£3 for the multichannel Kondo model and inCJ for the 
SU(N) x SU(f) Coqblin-Schrieffer model. We have gone 
further and have shown that other quantities, such as 
the entropy and in principle also the specific heat, can 
be calculated in the same way and for finite fields, as ex- 
plicitly described in Sec.IVa. However, more importantly 
we show how to obtain with the same numerical proce- 
dure uniformly accurate results for thermodynamics in 
the presence of an arbitrarily large level asymmetry, e/T, 
(or equivalently an arbitrarily large magnetic field, h/T, 
in the Kondo model). This was particularly important 
for obtaining the very low temperature thermodynamics 



of the asymmetric two-level system in this paper. Pre- 
vious techniques, for related TBA equations where the 
same technical difficulty arises, have either not been able 
to access large e/T (ft,/T)E3 or,L3, have required solving 
a separate set of approximate TBA equations suitable 
for large e/T (h/T). The latter requires matching the 
resulting thermodynamics at the boundary of the two 
temperature ranges. 

The iterative procedure to solve the TBA equations 
( |50| ) is as follows. An evenly spaced grid of 1000 points 
{A,, i = 1, . . . , 1000} is used in the interval -40 < A < 
+40, the functions £j(Xi),j = l,v are given initial val- 
ues Cj'(Aj) = (Aihand are then represented by least 
squares cubic splincsc3. The next iteration, £™ =2 (Ai), 
was obtained by evaluating the integrals to a relative 
error of less than 10 -5 at each point for each func- 
tion in turn. This procedure was repeated until the 

-N-l\ \ /lltiVll „• _ ^ 



, v with | \£j | 



vSi l£j(Ai)l 2 (when non-zero) reached machine preci- 
sion (approximately 10 -16 ). Convergence depended on u, 
with the procedure converging for a typical case after ap- 
proximately 30,44 and 61 iterations for a = (1 — 1/u) > 
1/2 and v = 3, 4 and 5 respectively. For fixed level asym- 
metry, e, we solved the above equations at each temper- 
ature. Similarly for a = 1/u < 1/2 the number of itera- 
tions required to reach convergence was typically 52, 79 
and 133 for v = 3, 4 and 5 respectively. 
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FIG. 17. Solutions to the TBA equations @ at 
aT/A r = 2~ 10 for a = 1/4 < 1/2 (a) and (b) (top) and 
a = 3/4 > 1/2 (c) and (d) (bottom) and for zero (left) and 
large (right) level asymmetry. The exponential decay in £i 
(£i,) for strong (weak) dissipation at large level asymmetries 
is clearly seen. 

A useful check on the correctness of the solution of the 
TBA equations was the value of the functions £j , Ej , Fj 
and Gj at A — * ±oo. These satisfy transcendental equa- 
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tions which can be solved analytically or numerically^!. 
We checked that our numerical solution for the above 
functions satisfied these boundary conditions at all tem- 
peratures and level asymmetries. 

In order to indicate the difficulty one encounters in 
implementing the above procedure for large asymmetries 
and low temperatures we show some results in Fig.]!?]. 
For zero asymmetry, the functions are well behaved and 
have variations on scales A ~ 1 so the interpolation 
scheme used above, which is crucial for carrying out the 
integrations accurately and ensuring the accuracy of the 
solution, is highly accurate. For large e/T we see however 
that for a > 1/2 (a < 1/2) the function £i (£„) decays 
exponentially to a vanishingly small value above (below) 
a characteristic rapidity A = B within an interval which 
vanishes as xo = gev /2T increasesc3. In this case a uni- 
form grid of points is inadequate and any interpolation 
will give wrong results. One could increase the density of 
points close to B, but a much simpler solution is not to 
solve explicitly for the functions £i (£„) but to solve in- 
stead for the quantities upon which they depend on (see 
Eq.©), namely S/ (A) = s * &(A) (s/(A) = s * &,_ 2 (A)) 
where s* is the integral operator defined in Eq. ( pl| ) . The 
latter functions s/(A) vary on a scale of O(l) and can 
safely be interpolated. From these one can then obtain 
£i (60 via its definition (|H^). The characteristic energy 
B is determined at each iteration and used to split up the 
integration ranges for those integrals containing £i (£„) 
so that the relative error of 10~ 5 for each integration can 
be maintained. 

The integral equations for the Ej , F.j and Gj (54 56|) , 
apart from being linear inhomogeneous integral equations 
as opposed to the non-linear TBA equations for the £j, 
have an identical structure to those of the latter. As 
above, for large xo and a > 1/2 (a < 1/2), one has to 
solve explicitly not for the E\,F\, G\ {E V ,F V , G v ), which 
vary rapidly near B as a result of the factors (1 — e - ^ 1 ) 
((1 — e~^)) appearing in the corresponding integral equa- 
tions (|54|~|5q), but for the quantities s * E2, s * F2, s * G2 
[s * E v -2-,s * F v ^2,s * GW2) with splitting of the in- 
tegration range for integrals involving the former. The 
same algorithm is therefore used as for the £j . Since this 
turned out to be a highly stable and accurate algorithm 
equally applicable without changes to both sets of equa- 
tion we did not consider the possibly simpler procedure 
of solving the linear equations for the Ej , Fj and Gj by 
Fourier transformation. This would involve numerically 
Fourier transforming the inhomogeneous terms, solving 
the resulting equations and then Fourier transforming the 
functions back to A-space. Computationally, this would 
be more efficient, but may not be as accurate as the pro- 
cedure we have used. In considering larger systems of 
coupled integral equations, v 3> 1, for which computa- 
tional time is a factor, the more efficient procedure may 
have to be implemented. 

The above technique, which resolves the problem of the 
exponential decay of certain functions to zero and thus 
allows the integrations at all points to be carried out 



with uniform accuracy, gives a controlled way of solv- 
ing the TBA equations. The thermodynamics is then 
calculated without taking any numerical derivatives. In- 
tegrating the functions £1 , E\ , F\ and G\ , weighted with 
,s(A + ///i) as described in Sec.IVa, gives the free energy, 
polarization, dielectric susceptibility and entropy respec- 
tively. Fig.|l8| illustrates how the above numerical diffi- 
culties were overcome for Gi(A, e, T) at low temperature 
and large level asymmetry. 



e/A=16, a=3/4 (v=4) 




FIG. 18. Gi(X,£,T) (solid line), and the integrand 
s(A + f/fi)Gi (A, e, T) (dashed line) of the expression for the 
entropy in Sec.IVa. The inset shows that the spike close to 
A = B ~ —2, due to the exponential decrease of £1 for large 
level asymmetry and low temperature, is a smooth structure 
on a small scale which is resolved by our numerical procedure. 
The integrated weight S(X) = J^ 4Q d\'s(\' + f/fi)Gi (A', e, T) 
(multiplied by 40000), giving the entropy, is also shown in the 
inset (long dashed line). The contribution coming from the 
spike is seen to be approximately the same as that coming 
from the broad feature. 



APPENDIX D: UNIVERSAL THERMODYNAMIC 
BETHE ANSATZ EQUATIONS FOR a < 1/2 

To derive the universal Bethe Ansatz equations one has 
to write Eq. jl9| ) in a form where the cutoff u> c , explicitly 
only present in the driving term Eq. (M) , can safely be 
removed. A simple analysis of Eq. (|19|) shows that for 
a > 1/2 in the A — > 00 limit all the e^'s tend to a finite 
value determined by T and e (or h) except for t\ which 
approaches e\ ~ — uj c and thus diverges in the lo c — > 00 
limit. This divergence is, however, harmless, since ei 
occurs only in the form ~ exp(ei(A)) at the r.h.s. of 
Eq. (p9[), and therefore all the 'interaction terms' behave 
nicely even in the tu c — > 00 limit. 

For a < 1/2 the situation drastically changes due to 
the sign change of the driving term. In this case one can 
easily convince himself that tj ~ oj c for j = 1, .., v — 1 
while t v ~ —lu c as A — * 00, i.e. the interaction terms 
blow up as one tries to remove the cutoff. The basic 
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idea to obtain the universal equations is therefore to 
transform Eqs. (|l9| ) into a form where only exp(— ej/T) 
(j = 1, ..,v — 1) and exp(e y /T) occur in the interaction 
terms. After a tedious calculation one finally arrives at 
the following equations: 



1 . 



ln(l + e^ /T ) - B jI * M 1 + er " llT ) 
Qj * ln(l + e e "/ T ) , (j = l,..,z/-l) 

u-l 

e u /T + Y,Qi * ln (l + e-^ /T ) 
X * ln(l + e e " /T ) , 



(Dl) 



(D2) 



where the integral operators K, Qj and Bji are defined 
by Eqs. 6.2.27 and 6.2.6 of Ref. [y| (Note that the last 
one of Eqs. 6.2.27 is written in Fourier space, and that 
the aj(A)'s occurring in the definition of Ajk are the func- 
tions given by the first line in Eq. 6.2.6 for all indices, 
j. Furthermore, in Eq. 6.2.6 the last Sjk term must be 
dropped in the definition of Ajk.) 

The driving terms, c/j in the equations above are de- 
fined as 



9u 



9] = l,..,v-2 



-9i>-\ 
2 



-w c artan e 



7rA/(^-l) 



tt r sinh Vj/ 

— + art an < 



■ sm 



' ttX \ 



2(^-1) 



(D3) 
(D4) 



where we corrected a factor of two with respect to 
Eqs. 6.2.28 of Ref. [ll|. Finally, approximating gj as 
gj = Cj exp{7rA/(z/ — 1)} and defining the dimcnsionlcss 
functions fj(X) = €j (A + ^-—^ In ^j-) /T one arrives at the 
following cutoff-independent universal equations: 



dj = ln(l + e^) - ^ S i; * ln(l + e 

2 = 1 

- Qj *ln(l + e^) , 
dv — + ^ * ln(l 

1=3 

+ K *\n(l + e v ") , 



(j = l,..^-l) 



(D5) 



(D6) 



where the 'universal' driving terms are defined as 

7TJ \ ( \tx 



<ij=i,..,i/-2(A) 
rf^ = — dj,_i = 



exp 



2(1/ -1) 
Avr 



exp 



f — 1 



i/ - 1 



These universal equations, apart from some sign changes, 
coincide with Eqs. 6.2.29 of Ref Finally, after some 
manipulations the impurity contribution to the free en- 
ergy can be written as 



T " f°° c%(A+^±ln^) 
F« = — — V / v< n Tk 

m.i. — ' 



d\ 



ln(f + e~ v ^ w )dX , (D7) 



where u, = + for j = 1 , . . , v — 1 and = -. The 
Kondo temperature Tk (renormalized tunneling ampli- 
tude) emerges naturally in course of the calculations, and 
is given by 

T K (a < 1/2) = 2De- vf / {v -^ 

= A r /a , (D8) 

which is identical with Eq. (|2^). 

We have checked that the solution of the above thermo- 
dynamic Bethe Ansatz equations gives identical results 
for the specific heat and static susceptibility as the much 
simpler equations from which they are derived. Fig. [l9| 
shows a comparison for the specific heat calculated using 
the two sets of equations. 



a= 1/3 



O 



0.40 



0.30 



0.20 



0.10 



0.00 
10 




FIG. 19. Specific heats, C(T)/T, for the symmetric 
two-state system (e = 0) at a = 1/3 from (a) the univer- 
sal form of the TBA equations (circles) and (b) the original 
untransformed TBA equations (solid line). 
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